Direct MD simulation of liquid-solid phase equilibria for three-component plasma 
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The neutron rich isotope ^^Ne may be a significant impurity in carbon and oxygen white dwarfs 
and could impact how the stars freeze. We perform molecular dynamics simulations to determine 
the influence of ^^Ne in carbon-oxygen-neon systems on hquid-solid phase equilibria. Both liquid 
and soHd phases are present simultaneously in our simulation volumes. We identify liquid, solid, and 
interface regions in our simulations using a bond angle metric. In general we find good agreement 
for the composition of liquid and solid phases between our MD simulations and the semi analytic 
model of Medin and Gumming. The trace presence of a third component, neon, does not appear 
to strongly impact the chemical separation found previously for two component carbon and oxygen 
systems. This suggests that small amounts of ^'^Ne may not qualitatively change how the material 
in white dwarf stars freezes. However, we do find systematically lower melting temperatures (higher 
r) in our MD simulations compared to the semi analytic model. This difference seems to grow with 
impurity parameter Qimp and suggests a problem with simple corrections to the linear mixing rule 
for the free energy of multicomponent solid mixtures that is used in the semi analytic model. 



I. INTRODUCTION 



The internal composition of White Dwarf (WD) stars has an impact on their evolution. Sedimentation of ^^Ne 
and the release of latent heat of fusion has been shown to delay cooling of WD stars by a few Gyr [T] . Such modest 
energy sources can have a large effect because the energy input from nuclear reactions is small. Therefore, the 
chemical energy input becomes relevant. Understanding these additional energy sources can allow for more accurate 
age determinations of stellar clusters. 

The interior of a WD is a Coulomb plasma of ions and a degenerate gas of electrons. As the star cools this plasma 
crystallizes. This crystallization has been observed in recent observations of both globular [5| and open star clusters 
[3] . The melting temperature observed in the 'Winget et al. results may also constrain the composition of the WD 
interior [3]. [Horowitz et al. only considered the contribution of carbon and oxygen. 

Much of the carbon, nitrogen, and oxygen, originally present in the star, is converted by nuclear reactions to the 
neutron-rich isotope ^^Ne. The larger mass to charge ratio of ^^Ne compared to ^^C and ^^O results in a release of 
gravitational energy as it sinks in the strong gravitational field of the star. This provides an additional source of 
energy during WD cooling :5,, 61. Sedimentation of ^^Ne can also affect the ^^Ni yield in type la SNe since the excess 
^^Ne in the core can modify the electron fraction [7]. 

Liquid-solid phase diagrams for multicomponent plasmas have been determined using Monte Carlo and density- 
functional techniques. Binary mixtures have been considered by many groups [8til2j. while three or more components 
have not been as extensively studied |13H17j . Often these works determine liquid-solid phase equilibria by comparing 
liquid and solid free energies that have been calculated separately. This approach may be sensitive to any small errors 
in the free energy difference between liquid and solid phases, while also providing no information on the dynamics of 
the phase transition. 
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Recently, we have performed direct two-phase molecular dynamics simulations of liquid-solid phase equilibria for 
carbon-oxygen mixtures in WD stars [4jll8). oxygen-selenium mixtures |18j . and for a complex 17 component mixture 
modeling the crust of an accreting neutron star [TS]. These simulations have both liquid and solid phases present 
simultaneously. This allows a direct determination of the melting temperature, and the composition of the liquid and 
solid phases from a single simulation. Systems with an arbitrary number of components can be modeled in this way. 

One must address potential systematic errors from finite size and non-equilibrium effects, however. Finite size 
effects can be important since an ion deep in the bulk liquid or solid may not be far from an interface. Larger system 
sizes must be used to address this. These two-phase simulations must also be run long enough to ensure that the 
phases come into thermodynamic equilibrium. This requires impurities to diffuse through the solid phase. However, 
diffusion in the solid phase is relatively fast since the ions have soft 1/r interactions, instead of hard cores, so that the 
ions can move past one another. We have extensively studied diffusion in Coulomb crystals in a recent paper |20| . 

If the systematic effects of a direct molecular dynamics simulation are addressed, then this method should yield 
accurate results. The systematic errors between the direct MD simulations and the free energy calculations are in 
principle very different so comparing the results directly should provide an important check on both methods. However, 
due to the lack of published results for free energy calculations of multicomponent systems, the only quantitative 
method currently available for comparing to these calculations is an extrapolation of their results for two-component 
systems [16]. 

In this paper, we perform MD simulations of C/O/Ne mixtures with 27648 ions for three different C/0 ratios as well 
as three different Ne concentrations for a total of nine systems. This choice allows us to see the neon dependence of 
the phase diagram across a range of C/0 ratios. We also compare our MD results to an analytic model extrapolating 
the results of two-component, free energy calculations. We discuss our MD formalism and analytic model in Section 
II, present results in Section III, and conclude in Section IV. 

II. FORMALISM 

We describe our molecular dynamics method in Section Ila, the algorithm to determine liquid vs solid vs interface 
in Section lib, and the theoretical model in Section lie. 

A. MD formalism 

The method used in these simulations is similar to our previous liquid-solid equilibria determinations [U 118) . We 
consider a three component mixture of ^^C, ^^O, and ^^Ne. The ions interact via screened Yukawa interactions 

Vtj{r) ^ —e 'I'', (1) 

T 

where Zi and Z^ are the respective charges of the two ions being considered, r is the separation between the ions, and 
A is the Fermi screening length, which for cold relativistic ions is 




Here, the electron Fermi momentum kp is kp = (Sn'^n^^ and a is the fine structure constant. The electron density 
Ue is equal to the ion charge density, rif, = {Z)n, where n is the ion density and (Z) is the average (by number of 
ions). Our simulations are classical and we have neglected the electron mass (extreme relativistic limit). Note that 
electrons can be non relativistic at lower densities and this will modify A. However, our results are not very sensitive 
to the exact value of A, see for example ref. [21j. 

One component plasma simulations can be characterized by the Coulomb parameter F = Z^e^ /aT. We characterize 
our multicomponent system using an average Coulomb parameter, 

F = (Z5/3)Fe, (3) 

where Fe — /afJT with the electron sphere radius — (3/47rne)^/'^. For our mixtures of carbon, oxygen, and neon, 

the values of ae and A are such that the dimensionless ratio n = a^j {Z)^^^ A ^ 0.4. Therefore, we expect the ground 
state to be a body-centered cubic (bcc) rather than a face-centered cubic (fee) crystal. S. Hamaguchi et al. pT finds 
K > 1.066, for a one-component plasma, in order for the system to be an fee crystal. 
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Time can be measured in our system in units of one over the plasma frequency ujp. For a one component plasma, 
the plasma frequency is 



M 



1/2 



(4) 



where M is the mass of the ion. For a mixture, we choose to define an average plasma frequency 
(47r(Z)2e2n/(M))^/^ 



B. Interface finding algorithm 



Determining whether a cluster of ions is a liquid or a bcc solid is simple when visually inspected, however this 
determination is difficult to obtain numerically. For an entire system, phase determination can be accomplished by 
computing the global order parameter Qq j^l] . In this work, we use the prescription laid out by ten Wolde et al. 
to determine whether individual ions are liquid-like or solid-like. 

For each ion i, an ion j is defined as a neighbor if it is within a given radius rmmi as defined by the first minimum 
in the pairwise correlation function g{r). The vectors fjj joining neighbors are called bonds. The direction of these 



vectors can be described by 0^ 
by a spherical vector qim{i), 



and 



in the frame of ion i. The local structure around ion i can be characterized 



(5) 



where N},{i) is the number of ions bonded with ion i and Yim{dij, 4'ij) is a spherical harmonic. 

These local order parameters are large in both the solid and the liquid. The global order parameter Qg is calculated 
from an average over all of the N ions. 



qsii 



1 ^ 



/6 = 



4 



-, 1/2 
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m— — 6 



(6) 
(7) 



This is large in the solid due to the fact that the qem(i) add up coherently. In the liquid, qemii) add incoherently, so 
Qq is near zero. This coherence is exploited to determine local order. For each qemii) a normalization is applied. 



qemii) 



qemii) 



_m— — 6 



1/2 



(8) 



A dot product can now be defined of the vectors qg for neighboring particles i and j. 



q6(«) ■ qeO') = E '?6™(*)%m(j)- 



(9) 



By construction, qg(i) ■ qg(i) — 1. 

We use the same criterion as ten Wolde et al. [53] for determining whether two particles are connected, namely 
qg(i) • qg(j) > 0.5. This criterion will be met for most of the bonds in the solid. In the liquid, two neighbors may 
be in phase and considered connected, but that is certainly not true for all of the neighbors. Therefore, we use a 
threshold on the number of connections to determine if an ion is solid-like or liquid-like. If an ion has more than 
seven connections, then it is tagged as solid-like. If it has seven or fewer connections it is liquid-like. For a perfect 
bcc crystal, the number of connections per ion is 14. 

Now that each ion is tagged as either solid-like or liquid-like, the interface in our two-phase simulations can be 
found. Deep in the solid, a vast majority of the ions within a certain radius of a given ion are identified as solid-like. In 
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FIG. 1: (Color on line) Automated interface finding algorithm results. Ions identified as being in the solid, liquid, or interface 
are colored black, purple, or pink, respectively. 

the bulk of the liquid a similar majority of ions are identified as liquid-like. Along the interface, there is a mixture of 
solid-like and liquid-like ions. For this reason, we tag an ion as being in the solid or liquid if a large majority (> 0.95) 
of ions within two lattice spacings are in the same phase. If this criterion is not met, then the ion is determined to be 
in the interface. 

Requiring a smaller majority, say > 75%, to define liquid and solid regions may lead to a thinner interface region. 
This could slightly increase finite size effects for the composition of the solid and liquid because we expect the interface 
region to have a composition, in general, intermediate between that of the solid and liquid phases. Therefore we choose 
an interface definition to yield a reasonably thick interface as shown in Fig. [l] Ions determined to be in the interface 
are found where one would expect them, along the border of the solid and liquid phases. Note that many of the ions 
determined to be part of the interface look to be part of the solid by eye. 



C. Semi analytic extrapolation of results from previous works 

There are a few calculations of three-component plasmas and in particular the carbon-oxygen-neon system [9l I13j , 
but it is difficult to make quantitative comparisons with these works. On the other hand, there are several quantitative 
results for one- and two-component plasmas [e.g. IHl HT]. Using the method described in Medin and Gumming |16) . 
hereafter MClO, we can extrapolate the one- and two-component results for comparison with the C/O/Ne system 
described in this paper. A full description of the extrapolation method is given in MClO; here we outline the basic 
algorithm. 

A multicomponent plasma (MCP) can be characterized by the charge Zi and fraction composition Xi — Ni/N of 
each ion species and the average Coulomb coupling parameter, F [Eq. ^^]. For a given Z = {Zi}, x = {xi}, and F, 
the equilibrium state of the mixture is fixed. This state, the state of lowest free energy, may be a pure liquid phase, 
a pure solid phase, or some combination of liquid and solid phases. For a given F, an MCP of composition x lies in a 
transition state between two phases of compositions a and b if it can be formed from a linear combination of the two 
phases and its free energy as a pure phase is greater than the total free energy of the combined phases; i.e., if 

Aa+{1-A)b^x (10) 

and 



Afa + (1 - A)fh < (11) 
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for some < ^ < 1, where fx is the free energy per particle of composition x, etc. The extrapolation method 
proceeds in two parts: First, the free energies of both the liquid and the solid phases of the MCP are calculated from 
the relevant one- and two-component values. Second, minimization is performed over (liquid and solid) composition 
space to find the mixture a, 6, A that gives the lowest total free energy for that composition, or to show that pure x 
is lower in energy than any other mixture with the same average composition. 

The free energy of the liquid phase of a multicomponent plasma is well described by the linear mixing rule (but see 
Appendix B of MClO); i.e., 

/;MCP(ri,xi,...,:r™_i) 



i=l 



/OCP(r,)+ln(.,A 



(12) 



where (Z) = '^"^i XiZ^ is the average ion charge and m is the number of ion species. The free energy of the solid 
phase of the MCP is 



/rcP(ri,xi,...,a;™_i) 



/f^P(r.)+ln(..A) 



i=l 

+ A/s(ri,a;i, . . . ,a;„_i) . (13) 



The expression we use for the deviation of the solid from linear mixing A/s is that given in Ref. 9J; another 
expression for Afs can be found in Ref. [TT] [see also [SJ US]. Here, fP'^^iTi) and /^'"^(Ti) are the free energies of a 
one-component plasma [as taken from, e.g. ITTl [T7) . 

For a multi-component plasma, assuming a is in the liquid state and b is in the solid state the minimization equations 
to solve are 

Ma) + ^{a) - a . V/K5) = fs{b) + §^{b) - 6 • V/,(5) , 

ie[l,m-l] (14) 

and 

Ma) - a • V/,(5) = fs{b) - 6 • V/,(&) (15) 



(see equations 28 and 29 of MClO). With Eqs. (14) and ( [T5| we must make a choice: Do we solve for F given the 
fraction < A < 1 of the mixture in state a, or do we solve for A given F? As will be discussed below, we choose the 
former for this paper. If we are given an average composition x and the fraction < ^ < 1 of the solution in state 
a (or the fraction 1 — A in state b), we can solve for F and the compositions of both the liquid and solid mixtures 
in equilibrium. We have 2m — 1 unknowns, ai, . . . , 6i, . . . , 6m- ii and F; but in addition to the m equations 



Eqs. ( 14 ) and ( 15 ) above we have the m — 1 equations 

Aai + {1 - A)bi ^ Xi , ie[l,TO-l]. (16) 

(note that ^ = 1, ^ 6^ = 1). 

Using the above method we generate results for the C/O/Ne system. These results, presented in Table |l| represent 
our best estimate of what Ogata et al. and other groups with similar free energy calculations would have found had 
they run their calculations for C/O/Ne, and in that sense allows us to compare our simulation results to those of 
previous works. For each result, the total composition {xi} and the liquid fraction of the system A are chosen to 
match values from a particular simulation run; given these values we calculate F and the liquid and solid compositions, 
Xq^Xqjx'^^ and XQ,XQ,x^g, respectively. Ideally, the value of A for each comparison is that of the liquid fraction 
used in the simulation run. However, the finite size of the simulation provides some ambiguity to the problem. In 
the simulations there is an interface region which is neither solid nor liquid that affects the total composition; in the 
semianalytic method discussed here the interface is assumed to be negligibly small (the system is assumed to be very 
large). To get around this ambiguity we effectively ignore the interface region and choose A to solve 

Ax^c + (1 - A)x'c = ■ (17) 

Note that a different value for A would be obtained had we chosen to solve for it in terms of xq or aiNc- However, 
A calcluated in terms of xq is very similar to that of xc and leads to differences in results of no more than 5%; A 
calculated in term of xng is also similar except when xnc ^ xc,o- 
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TABLE L Results using the theoretical model in Section [II C[ This approach finds the lowest energy liquid-solid mixture where 
Xi = Axi + (1 — A)xi. The total composition {xi} and the liquid fraction of the system A are chosen to match values from our 
simulation runs. 
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0.098 
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0.02 
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0.731 
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0.746 


0.227 


0.027 



Figure 2 of Ref. |4] helps show why we used the A values, and not the F values, from our simulations to make 
comparisons with previous works. Because of systematic differences between our direct MD simulations and the 
semianalytic extrapolation, for certain initial compositions and values of F (or equivalently, T/Tc in Horowitz 2010) 
the final liquid-solid ratios from these two calculations are qualitatively different. For example, there are regions of 
the phase diagrams at large F where mixtures are completely solid according to the extrapolation results but are 
liquid-solid according to the simulation results. By fixing the liquid fraction A instead of F we are guaranteed to 
find a liquid-solid mixture using the extrapolation method, though it will not necessarily be at the same F as the 
simulation result. 



III. RESULTS 

We now present results for liquid-solid equilibria for multiple mixtures of ^^C, ^^O, and ^^Ne. We intend to 
explore the effect of various trace neon concentrations xnc = 0.02, 0.10, 0.20 in different ratios of carbon to oxygen, 
3:1, 1:1, 1:3. For these compositions, the ratio of carbon to oxygen has a larger qualitative effect than different 
neon concentrations, so our runs will be grouped based on the carbon to oxygen ratio. 

A. Runs with 1 to 3 carbon to oxygen ratio 

Our first set of simulations had a carbon to oxygen ratio of 1/3 and wc used xnc = 0.02, 0.10, and 0.20. We 
randomly place 3456 ions in a simulation volume at a somewhat high temperature F ^ 150 relative to the expected 
melting temperature 175 < Fc ^ 300. We then make four copies of this system and this constitutes the initial liquid. 
We then take another 3456 ions at a lower temperature F ~ 300 and allow the system to crystallize. After evolving 
the solid system to allow it to roughly equilibrate, t ~ GOOOO/cjj,, we make four copies of the system to make the total 
solid initial configuration. The liquid and solid layers are then joined to make a total system size of 27648 ions. 

These initial configurations are not equilibrated for a number of reasons. First, since we combined a total of eight 
sets of initial conditions in building each large system, the energies along the boundaries may be high. There may be 
ions that are placed close together. The crystal is also not fully equilibrated with respect to structure and composition. 
These systems may take considerable time to fully equilibrate. 

We then evolve the system in small time steps. At ^ l/9uip. The temperature is adjusted to maintain roughly 
50% liquid and 50% solid. In Section lib we showed that the phase fraction is a difficult determination by eye when 
including the interface and generally have 50% liquid, 30% solid, and 20% interface ions. The system is then evolved 
for a total time of at least t — 2.5x 10^/cjp to allow the system to fully equilibrate and allow the impurities to diffuse 
through the system. These simulations were performed on the Cray XT5 system Kraken at the National Institute for 
Computational Sciences. 

Figure [2] shows the results of the first set of runs with xc/xq — 1/3, and x^ve = 0.20, 0.10, and 0.02. The number 
fraction of each ion species in the solid as a function of time is plotted in the left panel. The center panel contains the 
number fraction of the ion species in the liquid plotted as a function of time. The top right panel plots the fraction of 
the system that is solid. The interface region for all of our systems is roughly constant and about 20% of the volume. 
In order to determine liquid fraction one must subtract fs from ^ 0.8. We have also plotted F as a function of 
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time in the bottom right panel. The scales remain constant across Figures [2] [Sj and [4] in order to directly compare 
runs with different compositions. 

As the run starts, F was kept high (low temperature) in order to keep the badly non-equilibratcd solid frozen. As 
the system rapidly equilibrates, the temperature was then raised towards its final equilibrated value. The temperature 
then fluctuated at about the half percent level for the rest of the run. 




FIG. 2: (Color on line) (a) Composition of the solid region of an MD simulation vs simulation time tujp. The blue, red, and 
green lines indicates carbon, oxygen, and neon number fractions, respectively. The solid, dashed, and dotted lines correspond 
to simulations with xc/xo = 1/3 and xmb = 0.20, 0.10, and 0.02, respectively, (b) Composition of the liquid region of an MD 
simulation vs simulation time. These data are presented in the same fashion as the left-hand panel, (c) Fraction of the system 
that is solid fs vs tUp for xne = 0.20 (solid), 0.10 (dashed), and 0.02 (dotted). All of our simulations have roughly 20% of the 
simulation volume determined to be in an interface region. In order to compute the liquid fraction, on must subtract the solid 
fraction plotted from ~ 0.8. (d) T Ystujp for xnb = 0.20 (solid), 0.10 (dashed), and 0.02 (dotted). All of these simulations are 
done with 27648 ions. 

We expect the solid to be enriched in oxygen and the liquid to be enriched in carbon based on our recent carbon- 
oxygen phase diagram work [U I18j . The results shown in Figure [2] suggest that neon as an impurity does not 
qualitatively affect these results. 

B. Runs with 1 to 1 carbon to oxygen ratio 

Our second set of runs use a 1/1 carbon to oxygen ratio with xnc — 0.02, 0.10, and 0.20. We start these systems 
in a very similar way to the 3/1 carbon oxygen ratio systems. The history of these three runs is shown in Figure [s] 
We expect, once again, that the solid will be enriched in oxygen and the liquid enriched in carbon. As in the previous 
sets of runs, the addition of neon does not qualitatively change these results. 

C. Runs with 3 to 1 carbon to oxygen ratio 

Our final set of runs has an overall carbon to oxygen ratio of 3 to 1. We started these runs in a very similar way 
to all other runs described. The history of these three runs is shown in Figure |4] 
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FIG. 3: (Color on line) Compositions and phase fractions as per Fig. [2]except for initial compositions with xc /xo = 1/1- 



Our previous carbon-oxygen phase diagram suggest that there should be no difference in this ratio between the 
liquid and solid for carbon and oxygen. Now neon is enriched in the liquid phase in all cases. 



D. Collected results and supplemental run 



We now present the collected results from Sections III A[ III B III C The final compositions and final temperature 



of the runs are listed in Table |TT] These data are shown on a ternary plot in Figure [S] A point on this plot corresponds 
to the composition {xc,xo,XNe) subject to xc + xq + a^jve = 1- The number fraction of any given species is = 1 
at the labeled vertex and Xi = on the line opposite it. The solid symbols are the data from the MD runs and 
the open symbols are the model predictions. The fact that the squares and circles lie close to lines of constant neon 
concentration shows that the neon abundance does not strongly impact the shape of the CO phase diagram. 

The 3:1 carbon to oxygen ratio runs have the largest discrepancy between the MD runs and the theoretical model. 
In order to determine whether this is an actual difference between what was predicted and what was run or if the 
difference is due to statistical fluctuations, we performed another run for xc,o,Ne = 0.68,0.22,0.10 with twice as 
many ions, N = 55296, in a rectangular box with = ILj. = 2Ly, where are the box lengths in the three 
directions. We started this system with a biased composition as was done in Ref. [TH]- Instead of starting with the 
liquid and solid at xc,o,Ne = 0.68,0.22,0.10, we started with the composition predicted in the theoretical model of 
x'c ojy^ = 0.648,0.250,0.102 in the solid and x'-c^o,Ne = 0.712,0.190,0.098 in the liquid. The system was then run to 
t = 2.5 X 10^/a)p. These results are the blue triangles in Figurejs] Notice how even though the initial conditions agreed 
with the model predictions, this MD simulation still moved towards our previous N = 27648 ion MD compositions. 

We now test the agreement on the melting point, F, between our MD simulations and the theoretical model. Note 
that the model predictions were determined using our MD data for the liquid fraction, A. In general, the model 
predicts lower F values by 5-20%, see Table [ij than do our MD simulations, see Table In] and Fig. [e] This differ ence 
appears to be correlated with the impurity parameter Qimp, 

Q^rap= ^ X,{Z, - {Z)f . (18) 
i=C,0,Ne 
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FIG. 4: (Color on line) Compositions and phase fractions as per Fig. [2]except for initial compositions with xc /xo = 3/1. 



This parameter measures the dispersion in charge of the ions. For systems with a small number of impurities (small 
XNe) there is good agreement between model and simulations. However the difference is larger for systems with a 
large number of impurities (in general large Itvg)- This is illustrated in Figure [6] This difference suggests a possible 
problem with the correction to the linear mixing rule for the free energy of the solid phase that is used by the model, 
see A/s in Eq. (13 1 that is taken from Ref. [3]. For example the assumed bilinear form (see MClO equation 13) may 
be too simple. In future work we will calculate the free energy of solid C, O, Ne mixtures, using single phase MD 
simulations, and compare to the assumed linear mixing rule corrections. 



TABLE II: Composition of the total system, solid, and liquid phases for the nine MD simulations using molecular dynamics. 
The Coulomb parameter F is determined by the final temperature and the total composition of the run, see Eq. [3] 



xc 



Total 

Xo XNe 



0.20 0.60 0.20 

0.22 0.68 0.10 

0.24 0.74 0.02 

0.40 0.40 0.20 

0.45 0.45 0.10 

0.49 0.49 0.02 

0.60 0.20 0.20 

0.68 0.22 0.10 

0.68 0.22 0.10 

0.74 0.24 0.02 



ion 



27648 
27648 
27648 
27648 
27648 
27648 
27648 
27648 
55296 
27648 



232.75 
222.69 
211.02 
275.66 
248.93 
231.12 
307.58 
262.95 
261.01 
229.22 



Sohd 



''Ne 



0.125 0.661 0.214 

0.153 0.743 0.105 

0.170 0.810 0.0196 

0.320 0.460 0.220 

0.400 0.514 0.086 

0.424 0.556 0.0202 

0.585 0.224 0.191 

0.676 0.242 0.082 

0.652 0.252 0.096 

0.712 0.273 0.0149 



Liquid 



0.249 0.561 0.190 

0.271 0.633 0.096 

0.292 0.687 0.0201 

0.462 0.348 0.191 

0.481 0.409 0.110 

0.539 0.440 0.0204 

0.614 0.177 0.208 

0.684 0.200 0.116 

0.691 0.207 0.102 

0.758 0.218 0.0233 



IV. CONCLUSION 



We have performed molecular dynamics simulations to determine the influence of ^^Ne in carbon-oxygen-neon 
systems on liquid-solid phase equilibria. Both liquid and solid phases are present simultaneously in the simulation 
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Ne 



FIG. 5: (Color on line) Three component liquid-solid equilibria results. Vertices of triangle plot are 100% of the labeled 
ion species. For example, the vertex labeled C would correspond to a simulation with xq ~ i and the line connecting the 
other two vertices would be xc ~ 0. Small blue points are initial configurations. Filled circles (squares) show solid(liquid) 
equilibrated results. Open circles(squares) are from theoretical model. The blue triangles are results of the larger run where 
the upright(inverted) triangle is the solid(liquid). 
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FIG. 6: (Color on line) Comparison AF = Tmd — Fm between the MD simulation predictions for Coulomb parameter Tmd ~ F 
of Section HID and the model results Tm of Section II C versus impurity parameter Qimp, see Eq. (18 1. The carbon, oxygen, 
and neon number fractions, in %, of each point are labeled. 



volume. We identified liquid, solid, and interface regions in our simulations using a bond angle metric described in 
Section [II B[ This is the latest in a series of papers on liquid-solid equilibria both for complex multicomponent rp ash 
systems on accreting neutron stars |19| and for two component carbon-oxygen systems in cooling white dwarfs |31 118| . 

In general we find good agreement for the composition of liquid and solid phases between our MD simulations and 
the semi analytic model of Medin and Gumming [16 , see Fig. [5] The trace presence of a third component, neon, 
does not appear to strongly impact the chemical separation found previously for two component carbon and oxygen 
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systems. This suggests that the presence of small amounts of ^^Ne does not qualitatively change how the material in 
white dwarf stars freezes. This is consistent with assumptions in recent calculations of white dwarf evolution |24j but 
appears to be at odds with Segretain [T3]. 

However, we do find systematically lower melting temperatures by 5-20% (higher F) in our MD simulations compared 
to the semi analytic model, see Fig. |6] This difference seems to grow with impurity parameter Qimp. This suggests 
a problem with the simple corrections to the linear mixing rule for the free energy of multicomponent solid mixtures 
that is used in the semi analytic model. These simple corrections may have been fit to simulations involving only a 
small number of ions and where the distribution of impurities may not have been fully equilibrated. Alternatively, 
the model assumes corrections to linear mixing for the three component system can be written as a sum of pairwise 
terms and this assumption may be inadequate. 

To investigate these differences we are performing extensive single phase MD simulations of two and three component 
solid mixtures. These simulations involve 27648 or more ions and long simulation times. We are calculating radial 
distribution functions and static structure factors. In addition we are studying the diffusion of impurities and defects 
in these solid mixtures [20]. Not only do impurities occupy substitutional or interstitial positions in the crystal lattice, 
but the lattice itself can deform locally to accommodate impurities. This may be because the l/r Coulomb interaction 
has no intrinsic length scale. Therefore bond lengths are free to deform locally, depending on the charges of individual 
ions, as long as an average bond length is reproduced over large distances. In future work we will report on the 
structure of solid mixtures and on their free energies. 

This research was supported in part by DOE grants DE-FG0287ER40365 and DE-AC5206NA25396, by the National 
Science Foundation through XSEDE (eXtreme Science and Engineering Discovery Environment) resources provided 
by the National Institute for Computational Sciences under grant TG-AST100014, and by the Natural Sciences and 
Engineering Research Council of Canada and the Canadian Institute for Advanced Research. 
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